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We present and discuss a mathematical model for the operation of bilayer organic photo- 
voltaic devices. Our model couples drift-diffusion-recombination equations for the charge 
carriers (specifically, electrons and holes) with a reaction-diffusion equation for the exci- 
tons/polaron pairs and Poisson's equation for the self-consistent electrostatic potential. 
The material difference (i.e. the HOMO/LUMO gap) of the two organic substrates form- 
ing the bilayer device are included as a work-function potential. 

Firstly, we perform an asymptotic analysis of the scaled one-dimensional station- 
ary state system i) with focus on the dynamics on the interface and ii) with the goal 
of simplifying the bulk dynamics away for the interface. Secondly, we present a two- 
dimensional Hybrid Discontinuous Galerkin Finite Element numerical scheme which is 
very well suited to resolve i) the material changes ii) the resulting strong variation over 
the interface and iii) the necessary upwinding in the discretization of drift-diffusion equa- 
tions. Finally, we compare the numerical results with the approximating asymptotics. 

Keyuiords: Photovoltaics; drift-diffusion-roaction equations; finite clement methods; 
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asymptotic methods. 

AMS Subject Classification: 35K57, 35A35, 65N30 



1. Introduction 

The search for cheap, environmentally sustainable energy solutions has lead to in- 
tense interest and research in the area of organic photovoltaic (OPV) devices. Cur- 
rently, these devices have solar efficiencies which are significantly below those of 
modern inorganic semiconductor devices (i.e. at present 8.3% compared to 27.6% 
for "one-sun" Gallium Arsenide and 42.3% for concentrator celli^, currently lim- 
iting commercial implementation. 

A main difference of OPV devices (in contrast to inorganic PV) is the generation 
mechanism of free charge carriers, which typically occurs through the generation 
and dissociation of so called excitons, excited energy states created by incoming 
light. (We shall discuss excitons and the generation of free charge carriers in more 
detail below.) 

Many mathematical models have been proposed to study the behavior of bi- 
layer OPV cells as simple implementable organic devices: Refs. [H [Sj |3l [121 These 
models involve various approximations of the generation of free charges. A math- 
ematical basis for the these models comes from extensive literature on inorganic 
semiconductor models. For these devices, the standard macroscopic models couple 
two drift-diffusion-rcaction equations for the behavior of the free charge carrier den- 
sities (i.e. electrons and holes) in a sy stem with the Poisson equation governing the 
self-consistent electrical potcntial.^^ Much is known about such systems, including 
existence and uniqueness results for the steady-state (see e.g. Refs. HSl [14] and the 
references therein). The steady state is of particular interest for photovoltaic devices 
which arc expected to generate power on time scales much longer than the duration 
of the transient (i.e. the switching-on) dynamics. 

Drift-difFusion-reaction type models for OPV devices must crucially take into 
account the specific material properties of inorganic and organic semiconductor ma- 
terials, which involves in particular electric field dependent mobilities and specific 
recombination and dissociation rates (see Sec. [2| below). A second crucial amend- 
ment has to model the exciton dynamics and couple it appropriately to the original 
system. 

We shall thus consider the following evolution equations for the four main com- 
ponents, i.e. the concentrations of electrons (n), holes (p), and excitons {X) and 
the electric potential in the device (V), which must be calculated self-consistently. 

The evolution of the free charge carrier densities n and p shall be described by 
typical drift-diffusion-recombination equations and an additional source term due 
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to the excitons: 



dn 



V • (qD^Vn - qfinn\/{U + V)) - qRnp + qkdX 



(1.1) 




V • {qDp\/p + qHpp\/{U + V)) - qRnp + qkaX . 



(1.2) 



Here q denotes the positive fundamental charge, D„ , Dp and fj.p are the diffusion 
coefficients and niobihties of n and p, respectively, i?„p is the recombination rate 
of free electrons and free holes, and kd is the rate of dissociation of excitons into a 
pair consisting of a free electron and hole. Functional expressions for diffusion coef- 
ficients, mobilities and recombination/dissociation rates in organic semiconductor 
materials will be specified in Sec. [H 

For a bilaycr OPV device formed of two different organic semiconductor mate- 
rials (see Fig. [T|), the potential U (which is not an electric potential and doesn't 
contribute to the electric field) models an additional convection term that comes 
from the differences in the energy levels of the two materials. In an organic device, 
the energy levels involved in charge transfer are the Highest Occupied Molecu- 
lar Orbital (HOMO) and the Lowest Unoccupied Molecular Orbital (LUMO). The 
HOMO roughly corresponds to the valence band for classical semiconductors, and 
the LUMO resembles the conduction band. In general, the HOMO-LUMO gap cor- 
responds to the band-gap. 

For most organic materials, the HOMO-LUMO gap is too large for a photon 
to create free electrons and holes. Instead, it creates a bound electron/hole pair, 
a so-called exciton. In the bulk material, these excitons usually recombine without 
producing free carriers (with rate kr, see below). However, at an interface between 
two materials with suitable differences in the HOMO/LUMO properties, the exci- 
tons tend to align and split over the interface so that the electrons and holes are 
in separate (energetically favorable) materials. This effect is included in the model 
equations (jl.ip and (|1.2p by the term U, or more precisely, by the change in U over 
the interface region. 

Excitons are typically called polaron pairs (also referred to as charge-transfer 
states or coulomb bound pairs) when aligned across the interface and arc much 
more stable with respect to recombination. On the other hand, the aligned polaron 
pairs will cither be pushed together or pulled apart by an approximately parallel 
electric field. Accordingly the dissociation rate (kd) will be heavily field dependent. 
The combination of these two effects (reduced polaron pair recombination and field- 
driven dissociation) can lead to very high quantum efficiency (proportion of light 
creating free charges) for polaron pairs at an interface under an appropriately ap- 
plied external potential.^ Note that this is an internal quantum efficiency for polaron 
pair dissociation, not an external quantum efficiency, which further includes how 
many of the incoming photons actually create polaron pairs. 

Following the above discussion we shall assume that every exciton becomes 
immediately a polaron pair at the interface (see further explanation below). We 
therefore consider the following diffusion equation for the evolution of a combined 
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density of cxcitons and polaron pairs (to be identified by their position) , which lacks 
a convection term since the excitons are electrically neutral: 

dX 

— = V • (DxS/X) + cRnp + G-kdX- krX. (1.3) 

The coefficient G is the photo-generation rate for excitons, kr is the rate of geminate 
recombination of the excitons, and c represents the proportion of recombining free 
charges that form excitonic states (as opposed to recombination leading to emission 
of light, etc). 

Finally, we have the Poisson equation for the self-consistent electric potential: 

eoV-(e,Vy)-g(n-p-/ix/V,X), (1.4) 

where eg denotes the vacuum permittivity and is the relative permittivity of the 
material. Furthermore, represents the derivative normal to the direction of the 
material interface / and xi denotes the indicator function of the interface region 
while h is the separation length of the charge pairs of polaron pairs in the interface 
region. In fact, the extra term proportional to V^X is the field contribution due to 
the alignment of the polaron pairs at the material interface. The gradient results 
from taking the continuum limit of a sum of Dirac delta primes (corresponding to 
point dipoles). The exciton distribution X includes both polaron pairs and excitons, 
with their identity based on their location in the system. Polaron-pairs are split 
over the interface such that the electron and hole (although still bound) each lie in 
an energetically favorable material. This has the dual effect of aligning a sheet of 
dipoles (resulting in the aforementioned term in Eq. ()1^[) ) and of making polaron 
pair states extremely energetically favorable to excitons .1^^ We thus assume that any 
X in an interface region / are polaron pairs and that any X £ I"-^ (the complement 
of I) are excitons. In actual devices, organic interfaces are not sharp, and can be 
blended over a variety of length scales. For our model, we take / to be the region 
within d = Inm of the theoretical sharp interface. Note that if we take / to be a 
straight line parallel to the contacts of the device and assume homogeneity in the 
parallel direction, we exactly recover the model in Ref. [H 

We supplement the system (jl.ip to p.4p with the following boundary conditions. 
Dividing the boundary F into disjoint Dirichlet and Neumann parts F = F/j U F^r, 
we shall consider 

(V = Vo, n^nn, P^Pd on F,,, 

\v,l/ = 0, V,n = 0, V,39 = on Fw, " ^ ^ 

The boundary conditions for T ^ in Eq. (jl.Sp correspond to the conducting ends 
of the device (the left and right of Fig. [1]). The difference in the values for the 
potential V at the two Dirichlet contacts correspond to the potential difference in 
the device (in Volts) whose offset does not affect the dynamics of the device because 
only the electric field E = —\/V enters the other equations. The usual notation for 
a photovoltaic device follows the convention scheme of a diode, but in the working 
case current flows in the direction opposite to that of a usual diode. For this reason. 
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Fig. 1. Schematic of the 2-D device giving labels to the quantities used in the asymptotics and a 
graphical representation of the function U. For the numerics we take (7 to be a pieccwise quadratic 
connecting the values at either side so that U £ C^. Note that for the asymptotics the vertical 
direction of the first diagram is suppressed and the problem becomes 1-D. 



a consistent applied potential is given at the left boundary of the device. We thus 
take the potential at the right side of the device to be zero, and take the potential at 
the left side to be the potential difference. Note that this notation has the convenient 
consequence that the average field in the device has the same sign as the potential 
difference. 

For the boundary conditions of electrons and holes one has to take into account 
that the HOMO/LUMO levels of the device determine different energies for the 
electrons and holes implying that the concentrations of n and p at the metallic 
contacts can be very different. Thus, for majority carriers in their energetically 
favorable material {p to the left of the interface and n to the right as indicated in 
Fig. [T]), we take the boundary conditions given by Scott and Malliaras in Ref. ?IU[ 
For the minority carriers (n to the left and p to the right), homogeneous Dirichlet 
boundary conditions constitute a very good approximation as the Scott-Malliaras 
boundary conditions give values approximately six orders of magnitude smaller than 
typical concentrations (see Table 1). 

On the insulated boundary F^r we take homogeneous Neumann conditions for 
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all of the variables because no particles should be able to pass into the insulator and 
it should be electrically neutral. For the excitons, we take homogeneous Neumann 
conditions on the whole boundary F because the excitons are not be able to transmit 
from the organic material into the inorganic contacts. 

The results of this paper is as follows: In Sec. [2] we will scale the physical pa- 
rameters of the equations and detail the specific expressions for the mobilities, 
dissociation rates and recombination rates in an organic semiconductor material. 
Next, in Sec. 2] we perform an asymptotic analysis of the 1-D equations for some 
typical parameters. The asymptotic analysis quantifies approximately how the exci- 
ton production rate yields an electric current and gives an expression for the shunt 
resistance of the device. Moreover, we present the numerical system used for the 
2-D simulations in Sec. \5\ In Sec. |6] we plot the calculated concentrations in the 
device at three important points in the device operation - short circuit, optimal 
power point, and open circuit. We interpret these results in terms of their relation 
to the asymptotic results and furthermore show the current-voltage characteristic 
curves generated by both the numerical and asymptotic methods. 



2. Scaling and Models for Mobilities, Recombination, Dissociation 

We will scale (|l.ip - (|1.4p . The basic scaling introduces the following dimensionless 
quantities: 

V = = UtV, x = Lx, {n,p,X) = Nr{fi,p,X), ^ = MoM, (2.1) 

where is the thermal voltage, L is a characteristic length (usually on the order 
of the device length), Nr is a reference concentration, and fxo is a reference mo- 
bility. We assume the Einstein's relation D ~ ?7tM- Although this might not be 
justified for some organic materials, it greatly simplifies the equations and as of yet 
a generally accepted alternative model has not been found. We further introduce 
the dimensionless quantities: 

where T denotes a reference time scale. Note that the Debye length A is typically 
not a small parameter - see Table 1. 

For an OPV device under light we follow Ref. [3] and choose the Langevin re- 
combination rate = '?(^"+^p)"p which rescales with the reference time T as 
follows: 

_ 1 {fin + Mp)"P _ ~~ 
AT TT "P ~ \9 ~ CrUp. 

NrHoUr ^ fir- 
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Setting cj, = ccr this gives the foUowing system: 

dn 

'dt 
dp 

'dt 
dX 

'dt 

For the charge carrier mobihties and we shaU apply the Poole-Frenkel 
mode^S! postulating an exponential dependence on the square-root of the electric 
field E = -VV: 



n - 


~ hXi 










(2.3) 


V 


(/i„Vn - 




+ V)) 


— Crfip - 


F kdTX 


(2.4) 


V 


iUpWp + 


^ppV([/ - 




- Crfip + 


kdTX 


(2.5) 


V 




+ c'^.hp + 


GT- 


kdTX - 


KTX. 


(2.6) 



M„ = Mn(0)e^"Vli^l, ^p^^p(0)eWls|. (2.7) 

Moreover, the exciton dissociation rate kd shall be given by Ref. [3|as a function 
of M , which is a scaled square-root of electric field E: 

kd{M)-'-^{e^Hl^ir) + ij), 



where the superscripts represent positive and negative fields (with respect to 
the interface normal). 



Device Parameter values 

The physical scaling parameters we take are: 

V = UtV = .0258V F 
E = UtE/L = 2.58 X IG^V/mi? 

X = Lx = lOOnmi: 

n = Nrri = lO^^m"^ n. 

The values for the various constants are taken from Ref. [3] (except G which is 
converted to a spatial density instead of an area density). Table 1 collects all the 
used dimcnsionless parameters. 

The values of fix and kd^out are not easy to calculate physically and no consensus 
seems to exist in the literature for their values, but we take the given values based 
on reasonable estimates given in the sources listed above. 



3. Steady-State Equations 

For usual device operation, a solar cell will be producing current steadily for time- 
scales on the order of hours. Any transient behavior occurs over such short time- 
scales (i.e. microsecondjn]) that we neglect them for modeling the long-term be- 
havior and efficiency of the device. 
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Table 1. Values for physical parameters. 



(h/L) = .01 


(d/L) = .01 


er = 4 


A2 1.43 


T = .00386 


Esc = 13 


GT » 16990 


= 4 RJ .6987 


!7K) - Uixi) ^ 12 




ifxT 3.86 


fcr,oi.tT fa 3864 


fcd,in(0)TRi386 


fcd,in(+-Esc)T« 178 


ka,^„{-Esc)T ^ 2763 


M = 10-10 


(ixi/no) ~ -01 




/in(0) = 3 


7n = .788 


/^n(Ssc) ~ 53.3 


/.p(O) = 1 


7p = .153 


f^p{Esc) ~ 1.75 


n(xi),p(xo) as .04 


n{xo),p{xL) Si 4 X 10-'^ 





Note: Esc = 3.5 X 10'', corresponding to the short circuit potential difference of 
-.5V. The expressions hout and refer to the rate- values in- and outside of the 
interface region. 



Therefore, for the remainder of the paper, we will consider only the steady-state 
equations where wc drop all the tildes and absorb the time-scaling T into the rates: 

'A^V • {erVV) 

^ -V • (Ai„Vn - HnnV{U + V)) 
■ iUpVp + fipp\/{U + V)) 

We expect that this system of equations will have a unique steady-state solution. 
If we insist on the physically reasonable requirement that kd and fj, are smooth and 
bounded from above (corresponding to physical device saturation), then the first 
three equations fit very nearly into the standard semiconductor framework (see, 
for example, Ref. llSp . The only remaining difficulty is the exciton equation, but 
given n,p S D L°° as in the usual case, we see that for smooth /ix we have 
X £ W^'°° and the exciton terms in the first three equations should not pose any 
difficulty. Proofs of the existence and uniqueness and the corresponding results for 
the time-dependent are under current investigation and subject to an upcoming 
paper. 

Concerning the steady-state solutions, we emphasize that the light generation 
term G > implies a non-zero right-hand-side in the n and p equations of (|3.ip 
and thus non-constant drift-diffusion fluxes of n and p on the left-hand-side. Thus, 
we expect a steady flux of electrons and holes away from the interface with nearly 
negligible bimolecular recombination (due to the work-function considerations dis- 
cussed in the introduction). Even for G = 0, we will only recover constant stationary 
drift-diffusion fluxes in such particular cases such as fix = and Cr = j^J'^^ i which 
is not realistic given the considered parameters. As a consequence, the model sys- 
tem p.l[) is not designed to accurately predict the behavior of the device in the dark 
G = 0. 



n-p-^V,X 

— CrUp + kdX 

—Crnp + kdX 

c'^np + G — kdX — krX. 



(3.1) 
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In the following section, we shall investigate further the steady state solutions 
of ((XT|) for G > in a 1-D setting. 

4. 1-D Stationary State Asymptotics 
4.1. Large Applied Field Approximation 

We begin by considering the steady state equations in 1-D: 

-(Mr,"x - ^J'7ln{U + V))a; = -CrUp + kdX 
-{y-pPx + + = -CrUp + kdX 

-^ixX^x = c'^np + G - kdX - kr-X, 

where we have assumed that the and /ix arc spatially homogeneous. 

In reverse bias, which constitutes the operating state of an OPV bilayer device, 
a negative voltage Vdis < is given at the left boundary of the device, i.e. V{xo) = 
Vdiff , while zero voltage is given at the right boundary of the device, V{xl) = 0. Note 
that with this notation, Vdiff is the potential difference in the device. For working 
photovoltaic cells, this potential difference is the sum of the potential applied to 
the device and a built-in potential from the metallic contacts. 

With a constant relative permittivity we introduce the Debye length Xjj :~ 
X^Er ~ 5, and rescale the potential according to this potential difference, i.e. V — >■ 
I Vdiff I V which rewrites the Poisson equation from (j43| as 

V^oo = e (n- p- ^X^] , £ — TTITT — r V{xo) = -1, V{xl)=0- 

(4.2) 

Since in a working device a typical short circuit voltage is about —0.5 Volts,'^ which 
is |Vdift| ~ 19.3 in our units, we find that e sa 0.01. 

Hence (assuming for the moment that £j^Xx 1, where is the fraction of 
interface width to device width) we expect that the electric potential V is in zeroth 
e-order dominated by the potential difference and, thus that the electric field is 
constant in the zeroth order of e, i.e. 

V°ix) = Vdm-E'{x-xo), E':=-^^. (4.3) 

XL - Xo 

We shall see in the following that the approximation (|4.3p remains consistent 
with the assumption Sj^X^ <^ 1 after being carried over to the equations for the 
charge carriers and the excitons . 

However, first we quote from Ref. [20] the typical order of magnitudes of the 
boundary values of electron- and hole- densities in the working state of the device: 

n{xQ) < £, n{xL) = 0{e), p{xo) = 0(e), p{xl) < e. 

Next, we remark that with the zeroth order approximation also the mobilities 
and iJ,p are constant in zeroth. For the short circuit value {E'^l = 13, we calculate 
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that /x„ 50 and /Zj, 2 and thus l/(/z„/ip) w 0.01 ~ e. Thus, by rescaling 

/i„n — > n and jipp — > p, (4-4) 

we obtain with Cj. = 0(1) the foUowing rcscalcd charge carrier equations 

-{jix - n{y + U)a:)x = -0{e) np + kdX, hq := n{xo) < e, tll := n{xL) = 0(1), 

(4.5) 

-{Px + p{V + U):,)^ = -0{e) np + fcdX, po =: p(a;o) = 0(1), pl p(a;L) < £■ 

(4.6) 

Moreover, the rescaled equation for the excitons is 

-tixX,,^0{e)np + G~{kd + K)X, X,{0) = X,{xl) = 0, (4.7) 

where /xx ^ 1 is small, typically of order 0{e). 

Neglecting the quadratic 0{e)np term, we can solve this equation explicitly 
in terms of hyperbolic sines and cosines. This allows us to check the consistency 
the approximation (j4.3p and the rescaling (j4.4p with the underlying assumption 
£j;Xj; <^ 1: With kr.out ~ 4 • 10^ and fed, out ~ 1 outside the interface and with 
kr,in ~ 1 and kd,in ~ 3 • 10^ (for typical i?" < 0) inside the interface we have: 

\/ kd,in \/ kr.out J 

and thus e!}^\\X^\\L^^[,„^,^]) « 0(6^/4) if ^ = 0{e) and fix = 0(e). 
4.1.1. Zeroth Order Charge Carriers 

Now, we proceed to study the zeroth order solutions and p° of (|4.5p and (|4.6p . 
Introducing the zeroth order fluxes 

J,*] := - "V.., >/p := ~ipl + p\x). with = + [/, 

and neglecting the 0(e) recombination term, the equations (|4.5p and (|4.6p integrate 
immediately as 

^°(a;) = Jl^-F{x), J0(x) = + F{x), (4.9) 
F(x) := r kd{y)X'>{y)dy. 



Utilizing the Slotboom variables J° = e'^{n°e-'fi)^ and J^^ = -e"'^(p°e'^)a;, we 
can explicitly solve for nP and : 



n°{x) = noe^(^ 


)-Vo _|_ 


Jno f e'^^^'^-'^^y^dy 

Jxo 


= noe^^- 


)-Vo ^ 


p\x)=pQe^°- 


ip(a;) _ 


J Xq 


= Poe^°- 


ip{x) _ 
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where (po ■— fixo) and we define 



Jxo J Xq 



ip(y)-ip{x) 



Next, we can determine the parameters Jno^Jpo evaluating the boundary 
conditions u^xl) = n-L and p(xl) — Pl- Upon rearrangement, we have with ipL '■= 

^n[XL) ^ ^pKXl) 



and obtain exphcit formulas for 'nP and p*' : 

V '^n{XL) J "^nyXL) 



<^n{XL) 

(4.13) 

(4.14) 

which satisfy both boundary conditions since $„/p(x'o) = J'n/p{xo) — 0. 

These equations, although explicit, do not yet intuitively present the leading or- 
der contributions. However, using the fact that U is constant outside of the interface, 
one can explicitly calculate the integrals $„ and $p (see [Appendix A[ ). 

In the following, we use these formulas to identify the leading contributions to 
J^o and and thus and p". More precisely, we introduce the parameters 

1 gV(a:i)-¥'o^ }_ ._ ^•p{xr)-'p{xi) (4-15) 

S ' 77 

and observe that e'^^^'^'^''' = S^^ for the considered device geometry since we have 
2{xi~xo) ~ XL — Xr (see Fig[T|). With = —E^ + Ux and < as in the working 
device, these parameters are small: 5 « 10^"^ and 77 w 10^^ for = Egc with the 
numerical values given in Table 1. 

However, using these parameters directly in (|4.13p and (|4.14p does not directly 
yield insights into the behavior of and Because <^n/p{x) can change over 
many orders of magnitude, the behavior of and p° depends highly on competing 
exponential terms, none of which can be easily eliminated. One possible way to 
proceed splits the domain [xq, xl] = \xq, Xi] U {xi,Xr) U [xr, xl] into three areas: left 
of the interface, the interface, and right of the interface. For each of these areas, we 
would obtain formulas for and p^ of the form of (|4.10p and (|4.1ip in terms of 
the values •nP{xi), nP{xr) and p^{xi),p^(xr) and the boundary terms no, n^ and po, 
PL- Then, imposing continuity of n*^ and p'^ and continuity of the fluxes J^q and 
JpO at X = a;; and x = x,. would allow us to determine the values nP{xi), n'^{xr) and 
p°{xi),p°{xr). 
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However, a simpler way to proceed considers the zeroth order currents, which are 
sufficient to understand the produced electric current. Using the explicit formulas 
for $„ and in terms of 5 and ?/ given in [Appendix A| we have from ([4.12^ : 

JO ^ "'^-^+-^"(^^) . .4.16) 

^0 11/ i\ 1 \ f 1/ l^' 

W 1^5^ V ^ l) ^ EO-u^{e) 32 y- ^ 7j) ^ W ^ 5^} 

where Ux{0) denotes a mean- value of Ux within the interface. 

Because J„ + Jp = const, the sum J^q + J^q can be used to calculate the total 
current in the device (as predicted by the asymptotics). We shall plot and discuss 
the predicted relationship between the current and the applied field as asymptotic 
IV-curve in Sec. |6] (see Fig. Note that for the given plots we take U to be 
piecewise linear so that Ux{d) is explicitly deffiied. 



4.1.2. Short- Circuit Current 

As mentioned earlier, the sign of determines if the parameters 77, S are large 
or small. For the short circuit current, we have E'^ < and r],d ^ 0. In order 
to determine the lowest order terms of these equations, we must also calculate a 
more precise form of J-n/p{xL) and thus its order. At first we remark that F{x) = 
/^^^ kd{y)X^{y)dy = 0{lQ^) (since kd,inX^ = 0{G) = 0(10*) on the interface 
with width 10^^ and kd^outX^ = 0(10) outside the interface with device length 
xl^xq = 0(1)). Thus, one can verify with calculations similar to those below and in 
[Appendix A| that the exponentials of 'fi{x) will determine the leading contributions 
to the integrals Tnjpixh). More precisely, since ^(x) is strictly increasing in the 
< regime (recall that = —E^ + Vx) the leading order contributions derive 
from 

T^ixL) « / ' F{y)e^-~'^^y'>dy, J^p{xl) « / " F{y)e'^^y'>-'^-dy. 



Applying the definitions of Eq. (|4.15p . integration by parts and the mean value 
theorem yields 



r F{y)e^^-'^^yUy = ^ H F{y)e-^°^-'-y^dy 

Jxq V" J xo 

F{y) 



1 



1 F{xi) 1 F'{e,,) p -E^ix^^y) 

r]6^ EO TjS^ {E^Y 



Xo 
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for a mean value 0„ G (0,2;;). Together with similar calculations for J-p{xL) and 
9p G (xrjXi), we obtain 

v( ^^ 1 fcd.o.i^°(g») 1 (, ^\^n('H-^\ 

•^P(^i^ (^0)2 (1 - ) - + + ^('^ ) 

Taking only the dominating terms 0{rj^'^5^^) in (|4.16p and 0(1) in (|4.17p yields: 
Jl. « E^n, ^fL21^pM j;^ ^ k,.o^tXyp) ^ ^^^^^ 

« E^{no+PL) - %r(xO(0„) + X"(^?p)) - r k,X^{y)dy. 

The term E^{nQ+pL) gives a current which is proportional to the field, indicating 
that it represents a resistor. For a usual working device this term will be small 
(no + PL ~ 10~^ in Sec. [2]). Because the term does not depend on the length of 
the device, we interpret it as a shunt resistance (parallel to the device). We observe 
from Eqs. (|4.16p and (|4.17p that we have another current term with the form of a 
shunt resistance: rj5'^E^{nL +|Jo)- This factor is negligible in the short circuit case, 
but becomes large as we move to positive fields, where 5^ > l/rj. We investigate the 
total shunt resistance further in Sec. [6l 

The second two terms represent the interaction of the excitons on the system. 
For usual device parameters, kd^out ^ fcd.m (for the parameters given in Sec. [21 we 
have kd.out = 10~^fc(;,m)- Since X° = O(jr^) on the interface and X° = 0( ^ ^ ) 
outside the interface, we can neglect the contributions from outside the interface 
(recall that e (0, xi) and 6p e {x^, xi^)) compared to the contribution from within 
the interface. 

Thus the most important contribution to the current for negative fields is: 

^approx = - r ka,^nX\y)dy. (4.18) 

J xi 

This approximation works very well in the short circuit case (as examined in Sec.[6l 
in the discussion preceding Fig.|4]), where it replicates J^^j within 2% and is actually 
closer to the numerically calculated J. 

Note that J^pprox comes entirely from the J^q term. Because we have evaluated 
and Jp at the point xq, we expect that the hole current will be dominant. How- 
ever, using Eq. (|4.9p we can calculate and J^^, and that the only change from 
the currents at xo is that the —F{xl) term appears in the electron current. Thus, 
as expected, in the region where the electrons are favored, the primary contribution 
to the current comes from the electrons and the total current is conserved. 
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4.1.3. First Order Terms 

Since we have explicit forms for nP, p", and X'^, we can integrate twice to calculate 
(with an additional linear term to account for the two boundary values of V). 
Normally we don't plot V since it is generally dominated by its boundary terms. 
However, = —V^ can also be calculated explicitly (by a single integration) and 
we include the first-order field in the plots in Sec. \6\ The explicit form is not very 
illuminating and thus not written here. Note that we need to add a constant value 
to the field to insure that its integral gives the correct potential difference in the 
device (corresponding to the slope of the aforementioned linear term). 

In theory we could use the second order form for the potential to calculate 
the next order and solutions. In fact, this is more or less the essence of 
the iteration scheme outlined in Sec. [51 However, without a constant field, it is no 
longer possible to explicitly integrate the continuity equations (especially given the 
sub-exponential forms for the electron and hole mobilities). It would be possible 
to do these calculations numerically, but this would simply become a simplified 
version of our numerical iteration scheme, and thus we do not pursue this further 
in the asymptotic case. Instead, we will present another method for asymptotic 
calculations. 



4.2. Unipolar approximation 

The discussion of the previous section, in particular Eq. (|4.18|) . can be summarized 
by the statement that the electron- and hole- fluxes and are approximately 
constant outside the interface but feature a strong variation over the interface with 
a magnitude of F{xr) - F{xi) = kd,i„X°{y)dy. 

As a consequence, using the Eqs. (|4.13p and (|4.14p and similar calculations as 
in Sec. 14.1.21 and [Appendix A[ it follows that the zero-order electron- and hole- 
densities nP and vary also strongly over the interface. In fact, we obtain for the 
hole density p^: 

p\x,)^^^^^I^^^ + 0{5) + 0{p,) 

0, , _ F{xl)- F{Xr) , kd,^nX'>{e,^) , kd,outX'> {Oout) 



p'iXr) 



-f 0(77,52)+0(p^), 

where Oin S (a;;, x^) is an mean value within the interface, 6 out G {xr, xl) is a mean 
value outside of the interface and 77, 5 and p^ arc (as mentioned above) small in the 
working device. Then, since Ux{0) 3> 1 for 9 G [xi^x^) and because the dominant 
part of the integral F{x) comes from the interface region, we see that 

= 0(10^), in the working case E° < 0, 

P"[Xr) 

as a consequence of the work-function gap over the interface. A similar result can 
be obtained for n° giving n°(xr) ^ n'^{xi). 
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Note that in a physical device one might expect this ratio to be even larger. If 
we narrow the interface while keeping the work-function difference {U{xr) — U{xi) 
and the generation term given in ()4.18|) constant, we tend to increase the value 
of Ux(Pin)- Thus in some sense the smallness of p'^{xr) increases naturally from 
reducing the size of the interface. 

The previous estimates quantify the observation that a bilayer device is operated 
such that in one material, n I (here to the right of the interface) and in the other 
material p 1 (here to the right of the interface) because of the nature of the 
HOMO/LUMO levels of the different polymers. We can thus ask for an asymptotic 
simplification of system (|4.ip in the bulk-material away from the interface. 

In order to proceed with such an approximation, we observe furthermore that 

^'^'^o ~ 10 — 20, i.e. that the boundary value of the majority charger carrier 
p°{xi) given at the left of the interface dominates over the generation of holes from 
excitons in the bulk material left of the interface. In fact the above factor would 
even be larger in more realistic exciton models compared to ()4.1[) . which misses a 
term modeling a electrochemical trapping of excitons (due to stability of polaron 
pairs on the interfacJI|) yielding further decreased the exciton concentrations. We 
shall thus entirely neglect here the photogcneration term of charge carriers outside 
the interface {kd^outX^). 

Altogether, neglecting the minority charge carriers and excitons outside the in- 
terface, we consider the following simplified models for the majority charge carriers: 

-XlV^cc = XlE, =P, Vx-pE^ — ^, (4.19) 

in the p- material at the left of the interface (with the Debye length \\, = A^e^) and 

for the n-material at the right of the interface. Here Jpo and J„o are the constants 
arising from integration. 

For the rest of this section we shall focus only on the hole density p. The analysis 
of the electron equation is analogous. We remark that the Eqs. (|4.19[) constitute a 
generalization of the system considered in Appendix A of Ref. [5l where the zero- 
current case JpO = was analyzed. For the present case with non-zero current we 
assume first that /ij, is constant and exploit then a first-integral of the equation to 
obtain 

Ex-]:E^ = ^x + C, 

where Cp is another integration constant. Using the substitution x = —2u and 
E = Uu/y (see e.g. Ref. we obtain 



Vuu = {nu - '2Cp)y 
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where k = .f'''" . With the further substitution z = ^Tk (u — ) , this becomes 
the famihar Airy Equation, yzz = yz and we can write down the formula for y 
explicitly in terms of Airy Functions. Since E = yu/y, it is more convenient to 
write the form of ^ = — log {y) (derived independently in Ref. I16p : 



y = - loe 



X 2Cp 



2 



/3Bi 



X 2Cr, 



2 K 



(4.20) 



where we can express a and P in terms of the integration constants and the boundary 
values. Note that for the n material we take the substitution x — 2u and J„ positive 
and everything else proceeds identically. 

This asymptotic result gives an explicit form for the charge carriers in the bulk 
in 1-D. However, for bilayer systems away from the interface, the system generally 
approximates a 1-D system since the bulk material acts primarily as a charge- 
transport layer. Since we explicitly account for the self-consistent potential, this 
formulation works much better than our previous asymptotics when we have large 
carrier concentrations (especially for less negative potential differences). In theory 
we could combine these equations in the bulk with a numerical calculation of the 
interface, simplifying the computations. However, in practice, we can change the 
mesh (or grid) to account for the reduced complexity of the problem in the bulk 
(see Fig. (2) and the gain would not necessarily be significant. 

The main disadvantage of this formulation is that it relies upon the boundary 
data of the system. Since it doesn't model the interface, it is thus incapable of 
predicting the current in the device based on the device parameters, and requires 
J as a parameter. In fact, we need all of the boundary values poi -E'oj JpO in order 
to calculate the field and the hole concentration in the p-region of the device, but 
generally only po is given. The other values, and Jpo, both require calculating 
the solution of the whole problem or at least calculating a model of the interface and 
using asymptotic approximations of the boundary values Eq and Jpo ■ Alternatively, 
we can take the values p{xi), E{xi), Jpi and again calculate the values in the bulk. 
Either method gives extremely accurate results for the cases in which J <C 0, which 
we show in Fig. [4]- [6l 



5. Numerical Scheme 

Time- discretization 

We present a numerical scheme based on a Gummel-type iteration,l35i which is well- 
suited to deal with the nonlinear dependence of the equations for n, p, and X on 
the electric field [E ~ —WV). At each Gummel-iteration step we calculate first the 
electric potential V from the current state of n,p,X and then use this new electric 
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potential to update n,p,X. We ean write this as follows after the k step: 

- A^V • (e^V^fe+i) - Pfe - rife + ^VX^ 

- V • (Mri(i?fc+i)Vnfc+i - ^„(i?fc+i)nfc+iV(Vfe+i + U)) + CrUk+iPk = kd{Ek+i)Xk 

- V • {pp{Ek+i)Vpk+i + Mp(-Efc+i)pfc+iV(Vfe+i + U)) + CrUkPk+i = kd{Ek+i)Xk 

- —V • (/ijfVXfc+i) + {kd{Ek+i) + kr)Xk+i = c',.nkPk + G. 



Note that we use a semi-implicit discretization in time. Specifically, we implicitly 
discretize the flux and mass terms in (|5.ip . but use semi- implicit discretization for 
the recombination term CrUp and an explicit form for the dissociation term kdX. 
Thus n, p, and X can be effectively updated in parallel since they only require the 
previously updated values for V . Furthermore, each variable n, p, X is treated com- 
pletely implicitly within its equation, which greatly reduces numerical instabilities 
arising from the errors in the previous step. 

The iteration continues until a preset error between the two latest steps is 
achieved. In practical numerical simulation, we have found this Gummel iteration 
to converge well except for high positive values for the applied potential difference 
(which constitutes the non- working case). This case, which is characterized by op- 
posite convection terms coming from the work-function and the electric potential, 
leads to large values for n and p near the interface which prevent the Gummel 
iteration from converging. To bypass this difficulty, we add a damping parameter 
a which interpolates between old and new solutions. More precisely, denoting with 
Vj! , n'j, , p'j, , the solutions of the above Eqs. ()5.1|) . the damped scheme is: 



In particular, note that since V;!_|_i and Vk both satisfy linear boundary conditions, so 
does any convex combination of V^'_|_]^ and Vk - This is equally true for the boundary 
conditions for n,p, and X, which are fixed for a given potential Vk+i- 

Our numerical experiments in the case of positive applied potentials show that it 
is possible to optimize the convergence by tuning the damping parameter. In fact it 
turns out that using a small damping a ~ .01 for the first few (e.g. three) iteration 
steps serves to correct the initial guess and prevents the concentrations from going 
negative due to badly chosen initial data. Further optimization on the damping, 
such as a potentially dynamic choice of a, is a topic for further investigation. 

Space- discretization 

Our choice of a space-discretization has to take into account the two major diffi- 
culties inherent in our system. First of all, the material interfaces will cause (as 
discussed in the previous section) very strong changes in the densities n, p, and X 



(5.1) 



Vk+i=aVl+, + {l-a)Vk, 
Pk+i = ap'k+i + (1 - a)pk, 



Uk+i = a7^fc+l + (1 - a)nk, 
Xk+i = aX'f.j^^ + (1 - a)Xk. 



(5.2) 



February 6, 2012 1:35 WSPC/INSTRUCTION 

BrinkmanFellnerMarkowichWolfrani-OPV 



FILE 



18 

over and on the interface. Secondly, depending on the electric field, the equations 
can be either convection or diffusion dominated. Furthermore, we want our scheme 
to easily generalize to complicated 2-D interfaces. 

The Hybrid Discontinuous Galerkin (HDG) finite element method is well-suited 
to handle these challenges. The novel idea for the HDG method is to have sep- 
arate degrees of freedom on the elements u (an element in the usual sense) and 
on the facets (boundaries) of the elements up- We assume that our domain is de- 
composed into a mesh T consisting of triangles with a set of facets J-' consisting 
of the facets (edges) of these triangles. The space we use for each of the equa- 
tions is V = {{u,uf) ■ u e iJ^(fi) n H'^{T),up e L^(7')}. Our solution space is 
{V,Vf,n,nf,p,pf,X,Xf) ^ V = V x V' x V' x V' and we denote the approx- 
imate solutions to the four equations as UY,Un,Up,ux with corresponding facet 
functions uyj, Unf, Upf, uxf- The facet elements allow for upwind- type calculations 
for convection-dominated behavior, as well as allowing for large changes of behav- 
ior at material boundaries. The functions within the elements are approximated by 
polynomials of a given order. Both increasing the polynomial order and refining the 
mesh improve the properties of a solution after a given number of iterations at the 
price of larger computational costs. 

A solid introduction to HDG methods is Ref. [T31 specifically the first section 
on the general convection-diffusion case. The HDG method can more accurately 
be called a Hybridized symmetric interior penalty Discontinuous Galerkin method. 
The relationship to standard Discontinuous Galerkin methods runs quite deeply. 
The stabilization and symmetrization terms have close analogs in traditional DG 
methods, see for instance Refs. [71 |6l [H 

The specifics of the weak form of the equations are also given in Ref. [131 Here 
we shall recall only a few specific features. It is well-known that numerical methods 
for convection-dominated systems need to apply the correct upwinding. For our 
elements, this consists of choosing the facet elements for infiow boundaries, and the 
bulk elements for outfiow boundaries. This selection allows for calculation along the 
flow of the carriers in the usual upwinding sense. 

All the numerical examples consider the two-dimensional computational domain 
VL = [0, 1.5] X [0, 0.2] (in the scaled variables from SecllJ. In order to compare with 
the one-dimensional asymptotic analysis above, we take the domain interfaces (see 
Fig.[T]) to be straight vertical lines located at = 0, = 0.49, Xm = 0.5, x,. = 0.51, 
and XL = 1.5 respectively. 

The mesh generation was accomplished using the Netgen program created by 
Schoberl.fl^ In particular, the mesh generation automatically preserves the inter- 
nal interfaces allowing all coefficients to be defined piecewise depending on the 
properties of the local material. The numerics were done using NGSolve, a solver 
wrapper for the Netgen program. We also include d the inverse solver PARDISO to 
assist with the non-symmetric matrix calculations.^^^^I^See Ref. [22lfor more details 
about using the NGSolve program on convection-diffusion equations. 
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6. Discussion of Numerical Examples 

In the following wc present selected numerical examples of the steady-state device 
behavior of the organic photovoltaic bilayer device plotted in Fig. [T] In particular, 
we shall investigate the steady-states calculated from different applied potentials. 

Using the scheme outlined in Sec. Owe use the damped-Gummel iteration given 
in ()5.ip and (|5.2p with a = .01 for the first three steps (as a sort of preconditioning) 
and a = .6 until convergence. All of the plots show a mesh consisting of 720 elements, 
greatly refined near the interface, and use polynomial interpolants of order 7 for the 
bulk elements - see Fig. [21 Refining the mesh and increasing the interpolation order 
both give better accuracy, in particular for the steep concentration changes over 
the interface (see, i.e. Figs. SI [B] and [T]). Note that the geometry of the test cases 
is taken to be homogeneous in the y-direction and we plot only a sample x-cross 
section which can be compared with the one-dimensional asymptotics. A fully 2-D 
numerical example is given at the end of the section. 




Fig. 2. Mesh used for all of the plots in Sec.|6] 



6.1. IV Curve 

One of the most important characteristics of an OPV device is the effect of changing 
the applied voltage on the current in the device. We show this relationship by 
plotting the current as a function of potential difference by multiple simulations 
with changing boundary values. We pick values for Vdiff and then use these values 
to calculate the boundary values for n and p according to Ref. [20l Due to the built- 
in potential of the metallic contacts at the OPV device, we have that zero applied 
voltage leads to an approximate potential difference over the device of Vdis = —19.3. 

Fig. [3] plots the IV-curve for our typical OPV device (marked with +) with 
respect to the potential difference. The zero applied voltage point is approximately 
at the left edge of the graph. This referential IV-curve compares very closely with 
a second IV-c urve (marked with □) computed after replacing the stated boundary 
conditionJ^by homogeneous Dirichlet boundary conditions for each carrier at both 
Dirichlet boundaries. 
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Fig. 3. A comparison of tlic IV-curvcs obtained using tlic boundary conditions outlined above 
in Table 1 above and homogeneous Dirichlct conditions for the charge carriers and a plot of the 
generated power density for the usual boundary conditions. Note that the power corresponds to 
the areas of the rectangles with corners at the origin and the current at each applied potential. 



Both IV-curvcs show how the (absolute value of the) produced current decreases 
as the potential difference increases. In particular we observe a distinctive S-shape. 
This is an effect of the dependence of kd.in{E) on the self-consistent electric field 
E. The kink arises at the potential V « —3 for which the iJ-field touches zero 
from below at the interface. Note that the i?-field touching zero happens for a 
potential V smaller than zero because the holes on the left side of the device (i.e. 
the majority charges) increase self-consistently the negative applied £'-field towards 
the interface. This argumentation can be verified by plotting the IV-curvc obtained 
from a constant dissociation rate fed, in = kd.in{Esc)- A plot comparing this result 
to the asymptotic calculation is shown below in Fig. [TOl 

The power density generated by a device is given by P = ^(Mnt + Vdiff). As 
shown in Fig.|3l the power is zero for two specific characteristic points in the device: 
Vdiff — —Vint (short-circuit) and J = (open-circuit). A third characteristic point, 
the optimal power point, is defined by the voltage point with the maximal value 
of P. For potential differences outside of the regime shown in Fig. [3] the power 
generated is negative, indicating that the device requires power input to operate 
and is no longer in the working regime for a photovoltaic device. An upper bound 
for the maximum power is given by the product of the short circuit current and the 
open circuit voltage (these being the maximum current and maximum voltage for 
the device in the working regime) . The fill factor denotes the ratio of the maximum 
obtainable power and this upper bound.'^ Geometrically this is represented by the 
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largest possible rectangle contained under the IV curve and the smallest rectangle 
containing the IV curve. For the plots given in Fig. [3] the maximal power point 
occurs at a potential difference of F ~ —3 with a fill factor of 39.5%. This fill factor 
is in good agreement with other simplified organic devices.^ In general, fill factors 
vary heavily depending on the device, and are notably improved for multijunction 
devices. 



6.2. Electron- /Hole- and Exciton- Densities at SC, OPP and OV 

In this section, we study the details of the charge density and exciton distributions 
at particular characteristic points of an IV curve: i) Short Circuit (SC), ii) Optimal 
Power Point (OPP) and iii) Open Circuit Voltage (OV). The short circuit current 
is the current at zero applied voltage, or a potential difference of Vdifr = —19.3. The 
optimal power point is the potential difference for which the power is maximized, 
which is Vdiff = —3 in Fig. [3] The open circuit voltage is the point where J = 0, 
which occurs for Vdiff = 12. 

For each of these points, we will plot the numerical solutions for n (marked 
with □) and p (marked with x) in the same chart. On the same figure, we also 
plot the zero-order asymptotics derived in Sec. 14.11 Note that we have rescaled n 
and p by /x„ and in the asymptotics section, and must appropriately scale them 
back for the plots. Moreover, we show how the unipolar asymptotics derived in 
Sec. l4.2l reproduce the numerical solution in the bulk very well when complemented 
with the numerical data from the simulation of the interface (which is not part of 
the unipolar approximation). Specifically, we use the boundary data values p{xq), 
E(xo), and J for the holes and n{xL), E(xl) and J for the electrons. The unipolar 
approximation thus has the potential to reduce simulation costs for the dynamics 
of the bulk material. In particular we observe that in all three cases - SC, OPP 
and OV - the minority carriers (n to the left of the interface and p to the right of 
the interface) are strongly dominated by the majority carriers (as discussed above 
in Sec. 1121). 

Next, we plot the numerical solutions of the exciton concentration and the elec- 
tric field. The electric field is compared with the asymptotic approximation ex- 
plained in Sec. 14.1.31 The unipolar asymptotics do not include exciton behavior 
to compare with the numerical results. Furthermore, we refrain from plotting the 
unipolar approximation to the electric field outside the interface (which is discon- 
tinuous over the interface since the field is calculated independently on the two sides 
of the device). However, good agreement similar to the charge carrier densities can 
be achieved. 

We also list the numerically calculated current given at each voltage. The current 
predicted by the zero-order asymptotics is the sum of J^^^ and J^j given in Eq. (j4.16p 
and (|4.17p . However, because this predicted current ignores recombination, it only 
gives a good approximation in the short circuit case. 
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Short circuit 

The case with the largest (negative) potential difference Vdiff is the short circuit case 
(which has zero applied voltage). This is due to the built-in potential arising from 
the work functions of the anode and cathode metals. For the device parameters laid 
out in Sec. this occurs for Vdiff ~ —19.3. 

Since short circuit considers the largest negative potential difference it is thus 
the case for which the zero-order asymptotics derived in Sec. 14.11 are best suited. 
This is confirmed by Fig. 2] which compares the numerical solutions of n, p, X and 
E with the corresponding asymptotics. Note that we see a dramatic change in n, p, 
X at the interface, and that either electrons or holes are largely dominant on one 
side of the device. Moreover, note that the slope of the electric field is significantly 
larger on the left-side of the device due to the larger concentration of holes on the 
left-side compared to the smaller concentration of electrons on the right-side (Fig. 
2] plots n and p on different scales). 

The numerically calculated current (as listed in Fig.|4|) is J = —334.2. The cur- 
rent predicted by (|4.16|) . (|4.17|) is J° = —345.15. Although not in exact agreement 
with the numerically calculated current for the short circuit case, this constitutes a 
reasonable match for the zeroth order term. Furthermore, the lowest order current 
(in i5, rj) predicted by (|4.18|) is Japprox = —339.32, also in excellent agreement for the 
short circuit case. However, the current given by the asymptotics does not depend 
on the potential as dramatically as the simulations suggest, most likely indicating 
that the recombination term becomes increasingly important as we move away from 
short circuit. 

Another observed feature of Fig. |4] are the boundary layers for n to match the 
boundary data at xl and for p to match the boundary data at xq. The values of n 
and p at their majority side of the interface are predicted by the asymptotics to be 
n{xr) = 0.526 and p{xi) = 15.6, in good agreement with the numerically calculated 
values. We emphasize that Eq. (j4.18|) in Sec. 14.11 predicts that the concentrations 
of n and p at the interface are predominantly a consequence of the leading order 
fiux-changes over the interface ^appiox = ~ Xri"^ kd.inX^{y)dy and are essentially in- 
dependent of the boundary conditions po and n^. We can confirm this statement 
by using the interface values as the boundary conditions (i.e. n{xL) ~ 0.526 and 
p{xq) = 15.6). We plot the results in Fig. [5] without visible difference in the concen- 
trations of the charge carriers at the interface compared to Fig. ID In particular note 
that the current and the maximal values of n and p do not change significantly. 

Optimal Power Point 

A second key characteristic point for device operation is the optimal power point. 
This is the point where the power from the device is maximized. See Fig.[3]for a plot 
of this quantity. As shown above, this occurs for a voltage difference of Vdiff ~ —3. 
This is just below the voltage for which the resulting self-consistent field touches 
zero from below and the field-dependence of kd{E) causes a kink in the IV curve. 
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Holes (Numerics) 
Zeroth Order Asymptotics (p) 
Unipolar Asymptotic Fit (p) 
Electrons (Numerics) 
Zeroth Order Asymptotics (n) 
Unipolar Asymptotic Fit (n) 




0.6 0.8 
Position 



1.2 1.4 



Zeroth Order Asymproti 



7 



Fig. 4. Plots of n, p, X and E for the short circuit case Vdifi = —19.3. Note the different scales for 
n and p, and that the total variation of E is less than 10% of its value. J = —334.2. The electrons 
and holes are each dominant in one half of the device, and the asymptotic results provide good 
agreement to the numerical results in all cases. 



The concentrations n, p and X are plotted along with the E-tield in Fig. IB] In the 
plot of the £^-field we see a slightly more pronounced effect (compared to SC) of the 
excitons on the electric field (the small bump located at the interface), but again 
this effect is small compared with the overall changes in the field. 

Comparing Fig. |6] to Fig. |4] we see that n and p again are primarily located 
on their specific sides of the device {p to the left, n to the right) with an even 
larger change in concentrations across the interface. In contrast to Fig. 0] (which 
features boundary layers for n and p), the concentrations of n and p in Fig. [5] 
appear approximately linear in their majority sides. It is not clear whether this 
approximate linearity is a characteristic of the optimal power point or coincidental. 
Our numerical simulations, however, seem to suggest that this is indeed the case 
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16 I 1 1 1 1 1 1 1 1 0.6 




0.2 0.4 0.6 0.8 1 1.2 1.4 

Position 



Fig. 5. Plot of n and p for the short circuit case with artificial boundary conditions, = 0.526 
and po = 15.6. J = —332.528, compare these plots to those for the correct boundary conditions 
in Fig. U] The deviation of p from the constant solution results from the self-consistent potential. 
The effect is less pronounced for the electrons because the concentration is much lower. 

and we may conjecture that it is caused by a certain balance between drift and 
diffusion terms. Moreover, note that the scales for n and p remain different, but are 
considerably more balanced than in Fig. |4l 

Fig. [S] also plots the zeroth order asymptotics, which are still reasonable. More- 
over, we see excellent fits from the unipolar asymptotic behavior of the quasi-linear 
behavior of the n and p. The increasing discrepancy of the zeroth order asymptotics 
suggests the increasing importance of the nonlinearity of the self-consistent iJ-field. 
This can be confirmed by using instead of the numerical values E{xq) and E{xl) 
in the unipolar approximation which is then no longer nearly as accurate. 

Open Circuit 

The third point of importance for device operation is the open-circuit voltage. This 
is the point at which the current through the device is zero (analogous to an open- 
circuit in electronics). For our device this occurs at approximately Vdiff ~ 12. In 
physical units this corresponds to a potential difference of 0.31 Volts. Considering 
the internal voltage of our typical device, this corresponds to an open-circuit voltage 
of 0.81 Volts, larger than expected for a typical OPV. 

The plots in Fig. [7| show a strong alteration of the concentrations n and p 
compared to the previous Fig. 3] and [HI The hole concentration, p, changes from 
p{xi) ~ 370 to p{xr) ~ 1.4, over the interface, confirming that our numerical scheme 
is capable of capturing large changes across the interface (predicted to be 0(10'^) 
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n 0.2 0.4 0.6 0.8 I 1 .2 1 .4 0.2 11.4 n.6 O.S I 1.2 1.4 

Piisitiun Position 



Fig. 6. Plots of n, p, and E respectively for the optimal power point OPP Vdift = ~3. The resulting 
current is J = —241.585. The exciton contribution from Poisson's equation is visible as a bump in 
the field plot. We see that the zeroth order asymptotics begin to diverge from the numerical and 
unipolar solutions as the potential difference increases. 



by the zeroth order asymptoties). Note that for this case the maximum values of 
n and p are of about the same order albeit plotted on different scales (by a factor 
four). 

The behavior of the E-&e\d throughout the device is dominated by the large 
values of n and p with very little effect of the excitons at the interface. Note that we 
do not see a very good match of the asymptotic electric field, which is to be expected. 
In particular, since the zeroth asymptotics ignore the quadratic recombination term 
(crup), we are neglecting this loss term for n and p. This is particularly important 
within the interface, where the concentrations intersect at n sa p « 18 (recall 
that p is strongly decreasing from the left of the interface to the right while n 
is strongly increasing). This leads n and p to be overestimated. Note, however. 
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Fig. 7. Plots of n, p, and E respectively for the (nearly) Open-Circuit Vdiff = 12 case. J = 0.23. 
(The true Open Circuit Voltage is within 0.05 of 12, but further precision is not numerically 
relevant). Note the different scales for n and p. The qualitative behavior of the zeroth order 
asymptotics is nearly correct, but they no longer closely match the numerical results. In particular, 
the estimated value for p is off by a factor of 2 nearly everywhere as the majority carrier. 



that the asymptotics correctly reproduce the quahtative behavior of the carrier 
concentrations with exponential growth from the boundaries of the device to the 
interface. 

Finally we remark that the presented numerical scheme is capable of continuing 
further into forward-bias. However, the above model and in particular the param- 
eters are primarily focused on the working-case, and we do not necessarily expect 
good results for the strong forward-bias regime. In particular, as the number of 
polaron pairs on the interface increases (due both to a reduction in kd and an in- 
crease in cj,np), we observe that some of them diffuse away from the interface. This 
corresponds physically to polaron pairs becoming excitons, a transformation which 
is not generally physically observed. 
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6.3. Shunt-resistance and Asymptotics of the IV Curves 

A comparison of the IV curves in Fig.|3]indicates clearly that the nonlinear boundary 
conditions proposed in Ref. [201 appear to have very little effect on the IV-curve. The 
resulting difference is barely noticeable until we reach open-circuit after which it 
does grow as we pass into the forward bias regime. We can study this discrepancy 
more clearly by simulating the OPV device in the dark with G ~ 0. The resulting 
current is plotted in Fig. [S] and represents the shunt current in the device. Recalling 
the asymptotic currents Eq. (|4.16|) and (|4.17|) . we can see the increase in the shunt 
resistance for i? 3> (and thus (5 3> 1) proportional to the terms containing and 
Po- In the large S limit we recover a term of the form: (n^ +po)-E'- This has the 
usual form of a resistor following Ohm's law with (scaled) resistivity 1/(?t.l +Po)- 
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Fig. 8. An IV-curvc showing the similarity of the shunt current produced by injected carriers with 
no hght generation { X ) and the result of subtracting the IV curve for zero boundary conditions 
from the usual IV curve (+). See Fig. [3] for individual plots of these curves. The dotted line 
represents the zeroth order asymptotic fit for the current with G = 0. Note that the apparent 
noise for negative fields is a result of rounding errors in taking the diff'erence of the currents from 
two diff'erent simulations. 

The equivalence of this term to the difference between the nonlinear Dirich- 
let and homogeneous Dirichlet boundary conditions is heuristically clear from the 
curves plotted in Fig. [5] The difference shows exactly the effect of the boundary 
conditions given by Ref. dOl compared to zero B.C. For a device in the working 
case the boundary conditions are relatively unimportant, but as we move into the 
forward bias regime, these boundary conditions become more important. The given 
shunt resistance is reminiscent of the exponential growth of the dark-current for 
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a photovoltaic device in forward bias. Although both currents result from taking 
G = 0, our model is not designed to reproduce this dark current, which is normally 
generated by an equilibrium concentration term in the recombination rate - for 
instance Cr{np — riaoPoo)- With light generation, this term is negligible compared 
to the exciton contribution and thus is not included in our modeling assumptions. 
However, carefully choosing boundary conditions would allow our model to replicate 
the physically observed diode behavior in the forward bias regime. 

Note that the apparent noise for negative fields is a result of rounding errors 
in taking the difference of the currents from two different simulations. Choosing a 
smaller tolerance parameter for the convergence or calculating the current to higher 
precision (to eliminate 0(.l) errors for O(IOO) currents) will reduce the apparent 
noise. 

Next, Fig. HI plots the asymptotic IV curve predicted by the sum J„o + Jpo 
of Eqs. ()4.16|) and ()4.17|) . Although the asymptotic approximation is questionable 
for small potential differences (as stated previously) we do observe a qualitatively 
realistic IV-curve in Fig. [HI In particular we observe a similar open-circuit voltage 
and short circuit current as well as a significantly better fill- factor of 87%. The high 
fill factor for the asymptotic curves is expected, since the asymptotics neglect the 
recombination current, which should be the main cause of power loss for a bilayer 
device. Indeed, by ignoring the recombination term, the asymptotics describe what 
is called the photogenerated current, the idealized current generated by the device. 
This has interesting properties in and of itself, and this avenue of prediction is under 
further investigation. 



350 



Potential Diffcrcn 



Fig. 9. The IV Curve predicted by the asymptotics in Sec. 14.11 Note that it has a similar open 
circuit voltage and short circuit current, but has a better fill factor than the numerical simulations. 



A comparable increase in fill-factor is observed in Fig. [TUl where we plot a 
numerical IV curve assuming a constant dissociation rate fcd.m ^ fcr.in- Moreover, 
Fig. [TU] plots the corresponding asymptotic IV curve with good agreement. This 
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indicates that a primary loss of efficiency is due to the exciton recombination term, 
which becomes non-neghgible near the interface when approaching the open-circuit 
voltage. This quantifies the physical intuition that device efficiency is increased 
by forcing the polaron pairs to dissociate more efficiently into free carriers. For 
the numerical simulation the new fill factor is a much improved 60%, whereas the 
asymptotics are virtually unchanged from the standard model shown in Fig. O 
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Fig. 10. A comparison of the IV-curvos obtained using the numerical and asymptotic schemes with 
kd in = in{l^sc) = 2763 to be Constant. 



Constant k_d Numerics ■ 
Constant k_d Asymptotics 




We therefore plan to further study which parameter changes give the most im- 
provement in asymptotic device characteristics. By studying these various parame- 
ters, we shah aim to determine which devices are mostly likely to have the optimal 
characteristics, and then test whether they also imply an improved efficiency for the 
full numerical system. We expect to find interesting information in terms of device 
design and implementation on a fairly general level rather than optimizing our model 
to reproduce one particular pair of organic materials. Further work is underway to 
establish whether the predictions created by our model are directly applicable from 
a device design standpoint. The IV-plots shown here match qualitatively to those 
of realistic OPVs. A precise quantitative agreement will be challenging due to the 
large number of parameters in the problem, many of which are difficult to measure 
using standard physical techniques. 

Other further work is in progress addressing the rigorous mathematical theory 
of the system presented, as well as the possibility of deriving similar models for 
which the size of the interface tends to zero. 

Current interest also revolves around the development of intricate interface ge- 
ometries in solar cells. The most efficient solar cells do not use simple bilayer inter- 
faces, but generally use hetero junction interfaces which maximize the active area of 
the cell. Our model can be extended to the case of a bilayer with regular periodic 
interface, with the primary difficulty being the specific modeling assumptions for 
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non-straight interfaces. In Fig.[Tl]we show sample plots of the electrons and excitons 
for a non-uniform interface to demonstrate the capabilities of the numerics. 



l,556e-04 1.869e-.01 3.737e-01 -ff^giige-Ol 7,472e-01 




3,868e + 00 1.171e+01 l,956e + 01 ^^i^^Ql 3.525e+01 




Fig. 11. Sample numerical results for the electron and hole concentrations for a device with a 
non-uniform interface at short-circuit. The exciton concentration shows good confinement and 
demonstrates the shape of the interface. The strange peaks in n appear as a result of the given 
form of U and the concentrating effects of the electric field in the vertical direction. Additional 
modelling assumptions will be necessary to obtain physically realistic results. Note that away 
from the interface we see the same nearly-constant behavior with a boundary layer near the metal 
contact. 



Appendix A. Calculation of and ^p{x) 

First note that (px{x) = —E^ outside the interface and tpx{x) = —E^ + Ux{x) 
inside the interface. In addition tf(h) — ^(a) = ~E^{b~ a) outside the interface and 
(p{b) — f{a) = —E^{h— a) + U{h) — U (a) inside the interface. Note that the change in 
form of Lp forces different treatment of <p(y) depending on the region of integration 
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and a change in ip{x) depending on the value of x. We proceed by cases: 



Jxq " 

p ^-E\x-y) + U{x)-U{xi)^y _^ ^-E\x-y) + U{x)-U{y)^y 
J-^i ^-E\x-y) + U(xr)-U(x:)^yj^ 



W 



l-e 



-E"(x~Xn 



^°{x-x,) + U{x)~U(,Xi)^^(^^^^ 



^_g-B0(x-x,) + !7(x)-[/(x,) 



-B"-t/x(e) 



I _ ^-E\x-x^) 



Xo < X < Xl 

Xl '^C X Xf 

Xf X X ]^ 

Xq < X < Xl 
Xl X lCj^ 
Xr < X < Xl 



where we have included the U(xi) terms for clarity even though usually we take 
U{xi) = 0. In addition, the term Ux{0) for 9 € arises from using the mean 

value theorem after integrating by parts once. If we assume that U is piecewise 
linear, then = ig constant inside the interface and the integral is 



straightforward. In the same manner, we have: 



Xq < X < Xl 



^EOix-x,)-mx) + mx,)^^(^^^^ l^e- + Xl<X<X 



Xr < X < XL 



Recalling the definitions given in Eq. 14.151 r/ = e-^°(^''~^'^~'^(^'''+^''^'', S = 
gE (xi-xo) ^^Yl as 5^ — i^L-Xr) ^ Qa,n write the particular values $„(xl) 
and $p(xl) simply as 



%{xl) = ^ii5^{1-5) 



1 



1 / 1 



Note that the for < 0, we have (5^1, and the results are positive (as expected 
for integrating a positive quantity). 
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